Model  reduction  using  multiple  time  scales  in  stochastic  gene  regulatory  networks 

Slaven  Peles,*  Brian  Munsky,^  and  Mustafa  Khammasht 

Deapartment  of  Meehanieal  Engineering,  University  of  California 
Santa  Barbara,  CA  93106-5070 
(Dated:  August  28,  2006) 

Gene  network  dynamics  often  involves  processes  that  take  place  on  widely  differing  time  scales  - 
from  the  order  of  nanoseconds  to  the  order  of  several  days.  Multiple  time  scales  in  mathematical 
models  often  lead  to  serious  computational  difficulties,  such  as  numerical  stiffness  in  the  case  of 
differential  equations  or  excessively  redundant  Monte  Carlo  simulations  in  the  case  of  stochastic 
processes.  We  present  a  method  that  takes  advantage  of  multiple  time  scales  and  dramatically 
reduces  the  computational  time  for  a  broad  class  of  problems  arising  in  stochastic  gene  regulatory 
networks.  We  illustrate  the  efficiency  of  our  method  in  two  gene  network  examples,  which  describe 
two  substantially  different  biological  processes  -  cellular  heat  shock  response  and  expression  of  the 
pap  gene  in  Eseheriehia  eoli  bacteria. 


I.  INTRODUCTION 

Despite  the  great  progress  in  biosciences  in  recent 
years,  mathematical  modeling  of  intracellular  processes 
remains  in  its  infancy.  The  level  of  complexity  found  in¬ 
side  the  cell  is  such  that  intuition  obtained  from  other 
sciences  is  often  insufficient  to  develop  reliable  quantita¬ 
tive  models.  Many  cellular  processes  take  place  far  from 
equilibrium  and  on  timescales  longer  than  the  cell  repli¬ 
cation  cycle,  so  they  never  reach  the  asymptotic  state. 
Hence,  knowing  the  asymptotic  solution  alone  may  not 
be  sufficient  to  describe  the  system’s  dynamics.  Further¬ 
more,  characteristic  time-scales  in  intracellular  processes 
often  differ  by  several  orders  of  magnitude,  which  creates 
additional  computational  difficulties. 

At  the  essence  of  the  study  of  gene  regulatory  networks 
is  identifying  specific  regulatory  mechanisms  within  a 
complex  structure.  One  way  to  do  so  is  to  carry  out 
a  model  reduction  with  respect  to  the  process  of  interest. 
In  this  paper  we  propose  a  model  reduction  method  based 
on  singular  perturbation  theory  that  takes  an  advantage 
of  multiple  time  scales.  It  can  be  applied  in  the  study  of 
asymptotic  behavior  as  well  as  transient  processes.  The 
method  is  based  upon  different  mathematical  formalism^ 
and  in  many  aspects  is  superior  to  Monte  Carlo  based 
approaches^’^’^,  which  are  typically  used  to  treat  these 
problems.  The  accuracy  of  the  computation  is  known  a 
priori  and  can  be  adjusted  before  the  bulk  of  calculations 
is  carried  out. 

We  illustrate  our  method  using  two  examples  arising 
from  recent  experiments  with  Escherichia  coli  bacteria: 
we  analyze  the  pap  gene  regulatory  network  and  cellular 
heat  shock  response.  The  fact  that  the  two  processes 
occur  in  completely  different  biological  contexts  suggests 
that  our  method  is  suitable  to  study  a  broad  range  of 
problems. 

The  method  can  be  used  as  a  standalone  approach, 
but  its  applicability  is  significantly  extended  if  used 
in  conjunction  with  the  recently  proposed  Finite  State 
Projection  ^  Model  reduction  approaches  based  on  sin¬ 
gular  perturbation  theory  have  been  used  in  various  areas 
of  engineering  and  science^’^’^’^,  but  the  overwhelming 


complexity  of  some  biological  models  has  severely  lim¬ 
ited  applicability  of  these  approaches  in  some  areas  in 
biosciences.  The  Finite  State  Projection  retains  an  im¬ 
portant  subset  of  the  state  space  and  projects  the  re¬ 
maining  part  (which  can  be  infinite)  onto  a  single  state, 
while  keeping  the  approximation  error  strictly  within  pre¬ 
specified  accuracy.  The  resulting  finite  model  is  given  in 
an  analytical  form  and  allows  us  to  implement  further 
reduction  techniques,  such  as  singular  perturbation,  to 
a  broad  range  of  complex  systems.  As  it  reduces  sys¬ 
tems  to  manageable  sizes  at  little  or  no  cost,  the  Fi¬ 
nite  State  Projection  method,  combined  with  other  re¬ 
duction  approaches,  has  the  potential  to  fundamentally 
change  the  computational  approach  to  stochastic  biolog¬ 
ical  problems. 

This  paper  is  organized  as  follows:  in  Section  II  we 
give  a  brief  overview  of  some  computational  methods  that 
have  been  used  to  study  stochastic  gene  network  models. 
In  Section  III  we  describe  the  mathematical  details  of 
our  method.  In  Section  IV  we  provide  an  example  of 
how  our  method  can  be  applied  to  a  realistic  gene  net¬ 
work  problem,  and  in  Section  V  we  demonstrate  how  to 
use  time  scale  separation  together  with  the  Finite  State 
Projection  method.  Finally,  in  Section  VI,  we  discuss  our 
results  and  outline  prospects  for  further  research. 


II.  BACKGROUND 

Presently,  the  dominant  approach  to  modeling  gene 
regulatory  networks  is  to  describe  intracellular  pro¬ 
cesses  by  a  series  of  chemical  reactions  involving  pro¬ 
teins  and  RNA  molecules.  For  a  system  of  n  chem¬ 
ical  species,  the  state  of  the  system  inside  the  cell  is 
specified  by  copy  numbers  of  each  relevant  molecule 
X  =  (Vi,  V2, . . . ,  V^).  Often,  these  numbers  are  rela¬ 
tively  small  and  reactions  take  place  far  from  the  ther¬ 
modynamic  limit,  so  that  mesoscopic  effects,  most  no¬ 
tably  fiuctuations,  have  to  be  taken  into  account.  The 
state  space  of  the  system  is  not  continuous,  but  a  dis¬ 
crete  lattice,  where  each  node  corresponds  to  a  different 
X.  The  size  of  the  lattice  is  limited  by  the  maximum  pos- 
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sible  populations  of  the  n  chemical  species  in  the  finite 
volume  cell. 

At  the  mesoscopic  scale  one  describes  the  dynamics 
of  the  system  in  terms  of  the  probability  of  finding  the 
system  in  a  given  state  X,  rather  than  in  terms  of  tra¬ 
jectories  in  the  state  space.  The  dynamics  of  the  system 
can  be  modeled  by  the  master  equation  for  a  Markov  pro¬ 
cess  on  a  lattice^.  Although  respectable  attempts  have 
been  made  to  introduce  deterministic  mesoscopic  mod¬ 
els  for  chemical  reactions^^,  presently  stochastic  methods 
are  used  almost  exclusively  in  the  study  of  intracellular 
processes  at  the  mesoscopic  level. 

The  master  equation  describes  the  time  evolution  of 
the  probability  of  finding  the  system  in  a  particular  state 
X.  With  an  enumeration  X  ^  which  maps  each  pos¬ 
sible  state  to  a  single  index,  the  master  equation  can  be 
written  in  a  familiar  gain-loss  form^ 

KiPi  it)  -  WjiPi  (t)] ,  (1) 

where  pi  is  the  probability  of  finding  the  system  in 
state,  while  Wij  are  propensity  functions.  The  propen¬ 
sity  functions  define  the  probability  Wijdt  that  the  sys¬ 
tem  will  transition  from  the  to  the  state  during 
an  infinitesimal  time  interval  dt.  The  propensities  may 
be  obtained  from  the  chemical  reaction  rates,  which  of¬ 
ten  can  be  measured  experimentally.  The  first  term  on 
the  right  hand  side  of  the  master  equation  describes  an 
increase  in  the  probability  pi  due  to  transitions  to  the 
state  from  all  other  states  j,  while  the  second  term 
describes  a  decrease  in  pi  due  to  transitions  from  the 
state  to  other  states  j.  If  the  system  is  initially  found 
in  a  state  k,  the  initial  condition  for  the  chemical  master 
equation  can  be  written  as  Pi(0)  =  6ik,  where  6ik  is  the 
Kronecker  delta. 

The  solution  for  this  problem  is  the  probability  Pi{t) 
that  the  system  initially  found  in  state  k  will  be  in  state 
i  at  the  later  time  t.  If  we  define  Aij  =  Wij  —  6ij 
the  chemical  master  equation  can  be  written  in  a  more 
compact  form 

Piit)  =  Y^iiPjit)-  (2) 

j 

The  solution  to  the  chemical  master  equation  generally 
can  be  expressed  in  terms  of  evolution  operator  p(t)  = 
M(t,  O)p(O),  which  in  case  of  a  finite  A  can  be  written  as 

p(t)  =  exp(At)p(0).  (3) 

Solving  the  master  equation  at  first  seems  to  be  a  rather 
simple  problem,  as  there  are  many  efficient  methods  for 
solving  systems  of  linear  ordinary  differential  equations. 
However,  if  we  consider,  for  example,  a  process  involv¬ 
ing  three  proteins,  where  each  protein  comes  in  say  one 
thousand  copies  per  cell,  that  gives  us  up  to  a  billion  of 
different  states  and  a  myriad  of  possible  transitions  be¬ 
tween  them.  Carrying  out  calculations  for  such  system 


without  any  insight  about  its  biological  structure  would 
be  impractical  at  least. 

This  problem  may  be  ameliorated  by  using  a  Monte 
Carlo  type  of  computation^^.  The  idea  behind  this  ap¬ 
proach  is  to  start  from  some  initial  probability  distribu¬ 
tion  Pj{0)  =  Sjk^  then  choose  randomly  which  reaction 
will  take  place  next,  and  compute  the  new  state  j  where 
the  system  will  be  found  at  some  later  time  t.  The  prob¬ 
abilities  of  picking  a  particular  reaction  are  given  by  the 
propensity  functions.  The  hope  is  that  after  sufficiently 
many  calculations  like  this  the  histogram  containing  all 
outcomes  will  approximate  well  the  solution  of  the  chemi¬ 
cal  master  equation  p(t).  The  advantage  of  this  approach 
is  that  we  do  not  need  to  calculate  the  whole  matrix  A. 
Instead,  we  calculate  on  the  fiy  only  those  matrix  ele¬ 
ments  that  are  required  for  the  computation  at  hand. 
Furthermore,  this  method  is  broadly  applicable  as  it  re¬ 
quires  little  knowledge  about  the  details  of  the  system 
under  consideration.  It  has  been  demonstrated^^  that  in 
the  limit  case  of  infinitely  many  runs  the  Monte  Carlo 
solution  approaches  the  exact  solution  to  the  chemical 
master  equation.  Therefore  the  accuracy  of  the  compu¬ 
tation  can  be  increased  by  simply  generating  more  Monte 
Carlo  simulations. 

On  the  down  side  Monte  Carlo  methods  are  notorious 
for  their  slow  convergence^^,  and  the  amount  of  compu¬ 
tation  necessary  to  get  an  accurate  result  may  be  too 
large  to  be  completed  in  a  reasonable  amount  of  time. 
Also,  computers  cannot  produce  truly  random  numbers, 
so  one  has  to  generate  something  that  is  as  close  as  pos¬ 
sible.  Programs  called  random  number  generators'^  cre¬ 
ate  periodic  sequences  of  numbers  with  a  large  period, 
which  imitate  series  of  random  numbers.  If  the  period 
is  too  short,  periodic  patterns  will  create  numerical  arti¬ 
facts  in  the  calculation.  On  the  other  hand,  high  quality 
random  number  generators,  such  as  RANLUX^^,  take 
significantly  more  computer  processing  time  to  execute. 

Despite  their  shortcomings  Monte  Carlo  methods 
remain  an  important  tool  for  study  of  intracellu¬ 
lar  processes  and  variety  of  different  Monte  Carlo 
implementations^’^’^’^^’^^’^^’^^  have  been  successfully 
used  thus  far. 

An  alternative  approach  known  as  the  Finite  State 

Projection^d8,i9  been  proposed  recently  by  Munsky 

and  Khammash.  The  method  is  based  on  a  simple  obser¬ 
vation  that  some  states  are  more  likely  to  be  reached  in 
a  finite  time  than  are  others.  One  can  then  aggregate  all 
improbable  states  in  (2)  into  a  single  sink,  and  consider 
all  transitions  to  those  states  as  an  irreversible  loss.  This 
method  automatically  provides  a  guaranteed  error  bound 
that  may  be  made  as  small  as  desired  as  it  was  proved  in^ . 
With  some  intuition  about  the  system’s  dynamics,  such 
as  knowing  the  macroscopic  steady  state,  one  can  develop 
an  efficient  system  reduction  scheme.  It  has  been  demon¬ 
strated  for  a  number  of  biological  problems^ that  in 
this  way  the  system  (2)  can  be  reduced  to  a  surprisingly 
small  number  of  ordinary  differential  equations,  thereby 
dramatically  reducing  the  computation  time.  The  re- 
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duced  system  can  be  treated  analytically^  and  the  method 
does  not  require  computationally  expensive  random  num¬ 
ber  generation. 

By  discarding  unlikely  states  in  a  systematic  way,  the 
Finite  State  Projection  method  provides  for  a  bulk  sys¬ 
tem  reduction,  but  the  original  Finite  State  Projection 
stops  far  short  from  what  can  be  achieved.  For  example, 
the  method  does  not  consider  how  transitions  between 
the  remaining  states  take  place.  Transition  rates  be¬ 
tween  different  states  typically  vary  over  several  orders  of 
magnitude,  and  by  treating  them  equally  one  may  waste 
considerable  time  performing  computations  to  obtain  a 
precision  that  far  outstrips  the  model’s  accuracy. 

The  low  probability  transitions  occur  infrequently,  so 
the  processes  involving  them  will  take  place  over  long 
time  scales,  while  high  probability  transitions  correspond 
to  fast  intracellular  processes.  Different  time  scales  can 
pose  computational  problems,  as  the  system  of  ordinary 
differential  equations  (2)  becomes  stiff.  On  the  other 
hand,  depending  on  the  length  of  observation  time,  the 
system  can  be  further  simplified.  For  short  times  slow 
processes  may  be  neglected,  while  for  long  times  the  ef¬ 
fects  of  fast  processes  can  be  averaged. 

In  what  follows  we  introduce  a  computational  method 
that  addresses  these  shortcomings  by  taking  advantage 
of  multiple  time  scales  in  the  master  equation  to  sim¬ 
plify  the  system  of  equations  and  reduce  the  computation 
time.  This  method  is  in  a  sense  complementary  to  the 
Finite  State  Projection.  It  can  be  used  independently, 
but  significant  benefits  may  be  achieved  when  the  two 
methods  are  combined. 

III.  TIME  SCALE  SEPARATION 

In  order  to  define  a  proper  chemical  master  equation 
matrix  A  has  to  satisfy  some  general  properties.  Since  by 
definition  propensity  functions  Wij  ^  0,  all  off-diagonal 
elements  of  A  are  nonnegative.  For  the  same  reason,  all 
diagonal  elements  of  A  are  nonpositive. 

In  a  closed  system  the  probability  has  to  be  conserved, 
so  that  —  const,  for  all  times.  That  means 

^  X!  Pi  0  ^  ^  ^  (*)  =  0,  (4) 

i  i  j 

and  hence 

for  any  probability  distribution  p{t)  =  (pi(t), . . .  ,p7v(t)). 
Here  with  N  we  denote  the  number  of  all  possible  states 
where  the  system  can  be  found^^. 

Therefore  it  must  hold  ^  -  Aij  =  0,  i.e.  the  sum  of 
the  elements  in  each  column  of  A  must  be  zero.  In  other 
words  vector  1  =  (1,1,...,!)  is  a  left  eigenvector  of  A 
with  associated  eigenvalue  zero: 

I'^A  =  0.  (6) 


This  further  means  that  for  the  matrix  A  there  exists  at 
least  one  right  eigenvector  v  with  the  zero  eigenvalue, 

Av  =  0.  (7) 

The  eigenvector  v  represents  the  steady  state  probability 
distribution  for  the  system,  and  is  the  nontrivial  solution 
to  (2).  Furthermore  it  can  be  shown^  that  other  eigen¬ 
values  of  A  have  negative  real  parts  if  the  matrix  A  is 
irreducible,  i.e.  it  cannot  be  written  in  a  block  diagonal 
form. 

Note  that  we  are  interested  here  in  the  nontrivial  so¬ 
lution  to  (2),  which  exists  since  det  A  =  0.  There  also 
exists  a  trivial  solution  p  =  0,  but  we  can  discard  it  as 
nonphysical  since  it  does  not  satisfy  the  normalization 
condition  |p|  =  1. 

In  gene  networks  we  can  often  identify  clusters  of  states 
within  which  transitions  occur  quite  frequently,  while 
transitions  between  the  clusters  are  relatively  rare.  The 
chemical  master  equation  that  corresponds  to  such  a  situ¬ 
ation  has  a  nearly  block  diagonal  structure,  so  the  matrix 
A  in  (2)  can  be  written  in  the  form 

A  =  H^eV,  (8) 

where  is  a  block  diagonal  matrix  describing  transitions 
within  the  clusters,  matrix  V  describes  transitions  from 
one  cluster  to  another,  and  e  >  0  is  a  small  parameter. 

In  the  limit  case,  e  ^  0,  the  system  remains  in  one 
cluster  of  states  for  an  infinitely  long  time,  i.e.  the  prob¬ 
ability  for  the  system  to  be  found  in  one  of  the  states 
within  the  original  cluster  is  one.  Therefore,  same  as 
the  matrix  A,  each  block  of  H  should  define  a  proper 
master  equation.  Each  block  of  H  has  one  zero  eigen¬ 
value  with  associated  eigenvector  v^,  which  describes  the 
steady  state  probability  distribution  in  the  cluster, 
while  all  other  eigenvalues  of  the  block  have  negative  real 
parts. 

It  is  relatively  inexpensive  to  compute  the  full  eigensys- 
tems  for  the  smaller  blocks  of  H.  From  the  eigenvectors 
for  each  block  one  can  easily  construct  a  matrix  S  that 
diagonalizes  H: 

S-^HS  =  A,  A  =  diag(Ai,...,Ajv)-  (9) 

Matrix  S  has  the  same  block  diagonal  structure  as  H. 
This  procedure  is  further  simplified  if  some  blocks  of  H 
are  identical,  as  is  often  the  case.  We  label  eigenvectors 
and  eigenvalues  of  H  so  that  Re{Ai}  ^  Re{A2}  . . .  ^ 
RcIAat}.  The  first  m  eigenvalues  are  then  equal  to  zero 
{\^m  =0)  and  the  rest  have  negative  real  parts. 

In  order  to  keep  our  presentation  streamlined,  we  shall 
assume  that  for  all  negative  eigenvalues  |Re{Ai>^}|  ^ 
e.  This  is  always  satisfied  if  it  is  possible  to  make  a 
clear  distinction  between  fast  and  slow  reactions.  This 
assumption  can  be  relaxed  and  similar  results  obtained, 
as  we  shall  demonstrate  later. 

By  applying  now  S~^  to  both  sides  of  (2)  we  obtain 

X  =  +  eV^  X,  (10) 
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where  x  =  S'  ^p,  and  V  =  S  ^VS.  The  equation  above 
can  be  written  in  the  component  form  as 

N 

Xi  =  XiXi  +  e  VijXj.  (11) 

i=i 

From  singular  perturbation  theory  (see  Appendix)  there 
exists  a  near  identity  transformation 

x=(/  +  eG)y,  (12) 

that  removes  all  0(e)  terms,  which  depend  on  from 

the  first  m  equations  (i  ^  m).  In  other  words,  equations 
(11)  where  =  0  can  be  decoupled  from  the  rest  of  the 
system  by  a  coordinate  transformation  (12)  through  the 
order  0(e).  In  the  new  coordinates  the  first  m  equations 
become 

m 

Vi  =  e  y]  VijVj  +  0{e^).  (13) 

i=i 

By  truncating  O(e^)  terms  in  (13)  we  reduce  our  system 
of  equations  to  an  m-dimensional  problem.  The  reduced 
system  is  computationally  less  expensive  to  solve,  while 
it  still  approximates  well  the  dynamics  of  the  full  system. 
Because  (11)  has  a  stable  fixed  point  solution,  if  initially 
|x(0)  —  y(0)|  =  0(e),  then  for  all  times  t  >  0  it  holds 
W(.t)-y{t)\  =  0(e). 

Note  that  if  is  smaller  or  of  the  same  order  as  e,  the 
near-identity  transformation  (12)  and  its  inverse  intro¬ 
duce  corrections  to  the  equation  that  is  only  of  order 
O(e^).  Therefore  we  do  not  need  to  find  the  exact  form  of 
the  near-identity  transformation,  we  can  simply  truncate 
all  terms  containing  Xi^^n  from  the  system  (11). 

The  first  m  equations  can  be  solved  now  independently 
of  the  rest  of  the  system,  and  their  solution  can  be  written 
as 


m 


where  y'  is  m  x  m  principal  submatrix  of  V,  with  ele¬ 
ments  Vij^rn-  In  many  cases  of  interest,  solving  (13)  is  a 
manageable  problem,  unlike  getting  general  solution  for 
the  chemical  master  equation  (2).  Since  in  the  long  time 
limit 

lim  Xi^rn(t)  =  0(e),  (15) 

t^OO 

as  we  show  in  the  Appendix,  we  claim  that  from  the 
solution  to  the  truncated  system  (14),  we  can  easily  con¬ 
struct  an  approximation  to  the  full  solution  of  the  chem¬ 
ical  master  equation  (3).  To  do  so,  we  first  define  an 
evolution  operator  V(t)  such  that  Vij{t)  =  [ex.p{eV't)]ij 
for  i,jf  ^  m,  and  Vij{t)  =  0  otherwise.  In  a  block  matrix 
form  this  is 

i>(t) = ( “V s )  ■ 


The  price  we  pay  for  simplicity  here  is  that  V(0)  is  not  an 
identity  matrix,  so  the  initial  condition  ^i>m(0)  also  gets 
truncated,  and  that  truncation  error  is  generally  larger 
than  0(e).  The  transient  time  T(e),  after  which  the  solu¬ 
tion  of  the  truncated  system  is  0(e)  close  to  the  solution 
of  the  full  system,  can  be  estimated  from  the  leading 
nonzero  eigenvalue  as 

T(e)^lne/A^+i.  (17) 

If  the  time  scale  separation  in  (8)  is  done  accurately,  this 
transient  is  negligible  for  all  practical  purposes.  However, 
time  scale  separation  in  a  large  system  is  not  always  ob¬ 
vious,  and  may  be  error  prone.  We  shall  further  discuss 
this  issue  in  Section  III  A. 

Finally,  by  performing  the  inverse  S  transformation  on 
V(t),  we  obtain 

V(t)  =  SV(t)S-\  (18) 

which  leads  to  the  0(e)  approximation  to  the  solution  of 
the  chemical  master  equation  (2),  that  is 

|p(t)  -  V(t)p(0)|  =  0(e),  (19) 

for  all  times  t  >  T(e). 

Note  that  neither  V  nor  V'  are  generators  for  the  evo¬ 
lution  operator  V,  so  their  eigenvectors  cannot  be  used 
directly  to  compute  the  steady  state  probability  distri¬ 
bution  for  p(t). 

A.  Computational  algorithm 

One  should  note  that  due  to  the  truncation  of  V  only 
contributions  of  the  first  m  columns  of  S  and  m  rows 
of  S~^  affect  the  approximate  solution.  As  a  result  the 
computation  can  be  greatly  simplified  -  instead  of  cal¬ 
culating  full  eigensystems  for  each  block  Hi,  it  suffices 
to  find  only  the  eigenvectors  associated  with  the  zero 
eigenvalues.  Instead  of  S  we  use  the  N  x  m  matrix 
whose  columns  are  made  up  of  the  right  eigenvectors  of 
H,  while  instead  of  S~^  we  use  the  m  x  N  matrix 
whose  rows  are  made  up  of  the  left  eigenvectors  of  H. 
Note  that  the  left  eigenvectors  are  always  if,  provided 
all  |v^|  =  1,  so  the  matrix  is  obtained  at  no  compu¬ 
tational  cost.  The  accuracy  of  the  calculation  is  known 
a  priori  to  be  0(e)  for  all  t  >  T(e). 

To  improve  the  reliability  and  robustness  of  our  calcu¬ 
lation,  we  can  optionally  add  a  transient  time  check  to 
our  algorithm.  To  do  so  we  first  need  to  find  eigenval¬ 
ues  for  all  blocks  Hi.  This  comes  at  a  relatively  small 
computational  cost,  and  can  be  performed  before  all  the 
other  calculations.  The  transient  time  needed  to  obtain 
the  desired  accuracy  is  estimated  from  the  leading  nega¬ 
tive  eigenvalue  A^+i  according  to  (17).  If  the  transient 
is  too  long,  that  can  be  remedied  by  expanding  matri¬ 
ces  and  to  include  the  right  and  left  eigenvec¬ 
tors  corresponding  to  A^+i,  respectively.  The  transient 
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time  is  then  governed  by  next  negative  eigenvalue  A^+2- 
This  procedure  can  be  repeated  until  the  desired  accu¬ 
racy  is  achieved,  thereby  sacrificing  computational  time 
for  precision.  Note  that  in  this  case  the  right  eigenvectors 
corresponding  to  nonzero  eigenvalues  cannot  be  obtained 
trivially. 

By  performing  this  test,  we  also  ensure  that  condition 
|Ai>m|  ^  e  is  satisfied.  Eigenvalues  of  H  that  are  0{e)  or 
smaller  will  result  in  long  transient  times.  By  expanding 
transformation  S  to  include  eigenvectors  corresponding 
to  these  eigenvalues,  we  essentially  treat  them  as  if  they 
were  part  of  V.  This  procedure  adds  robustness  to  the 
method  with  respect  to  separating  fast  and  slow  reactions 
in  (8). 

Let  us  summarize  how  to  implement  our  method: 

1.  Specify  problem  parameters.  If  necessary  apply  a 
finite  projection  to  the  full  state  space.  Use  propen¬ 
sity  functions  and/or  physical  intuition  to  separate 
H  and  V. 

2.  Find  the  eigenvalues  of  the  uncoupled  system,  and 
identify  “slow”  ones  with  respect  to  a  preset  tran¬ 
sient  time  T(e). 

3.  Find  right  and  left  eigenvectors  corresponding  to 
the  slow  eigenvalues  and  construct  rectangular  ma¬ 
trices  and  S^. 

4.  Calculate  k  x  k  matrix  V'  =  where  k  is 

the  number  of  slow  eigenvalues. 

5.  Compute  k  x  k  matrix  exp(eUT). 

6.  Perform  the  inverse  transformation 

S^expi^^V't)  =Vit) 

in  order  to  obtain  the  approximation  to  exj[){At)  for 
all  times  t  >  T(e). 


B.  Example 


Let  us  illustrate  this  technique  with  a  simple  example. 
We  assume  two  weakly  interacting  systems  that  can  be 
found  in  three  different  states  each.  We  choose  matrices 
H  and  V  in  an  arbitrary  way,  with  the  only  constraint 
that  they  define  a  master  equation.  In  our  example 


H  = 


Hi  0 

0  H2 


is  a  block  diagonal  matrix  with  blocks 
/  -4  2  4  \ 

Hi=  \  1-2  0  1  and  iLs  = 


-6 


(20) 


2  -3 


0  -2 


(21) 

We  find  that  blocks  Hi  and  H2  have  one  zero  eigenvalue 
each,  with  corresponding  right  eigenvectors  vi  =  (4,  2,  3) 


FIG.  1:  Comparison  of  the  approximate  and  the  exact  so¬ 
lution  to  the  master  equation  in  Section  IIIB.  The  initial 
probability  distribution  is  Pi{0)  =  621.  The  transient  time  is 
estimated  to  be  T(e)  =  Ine/As  =  1.96  for  e  =  0.01,  and  is 
denoted  by  the  vertical  line  on  the  graph. 


FIG.  2:  1-norm  error  in  probability  distribution  for  the  trun¬ 
cated  solution  versus  e.  For  each  value  of  e  we  have  randomly 
generated  50  matrices  H  and  U,  so  that  every  H  -\-eV  defines 
a  proper  master  equation.  Each  matrix  H  has  between  2  and 
6  blocks  and  each  block  has  size  between  2  and  21.  The  el¬ 
ements  of  H  and  V  are  randomly  generated  from  a  uniform 
distribution  between  0  and  1.  The  probability  distributions 
were  calculated  after  time  t  =  2T(e)  =  21ne/Am+i- 
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and  V2  =  (3,2,6).  From  these  eigenvectors,  we  assemble 
the  matrix  5'^, 


/  4/9  0  \ 

2/9  0 

3/9  0 

0  3/11 

0  2/11 
0  6/11/ 


(22) 


The  matrix  composed  of  left  eigenvectors  of  Hi  and  H2 
is  similarly  used  to  form  5'^, 


/  1  1  1  0  0  0  \ 

0  0  0  1  1  1  y  • 


(23) 


In  our  example  the  coupling  matrix  is 

/-8  0  0  5  32\ 

0  -5  0  2  3  1 

0  0  -12  4  6  2 

^  “  4  2  3  -11  0  0 

1  2  5  0  -12  0 

V  3  1  4  0  0  -5  / 


(24) 


To  get  the  equations  for  the  slowly  changing  variables 
(13),  we  calculate 


V'  =  s^vs^ 


-87/11  78/11  \ 
29/3  -26/3  )  ■ 


(25) 


Next,  we  calculate  the  evolution  operator  for  the  trun¬ 
cated  system,  ex.p{eV't),  and  perform  the  inverse  S- 
transformation  to  obtain 


V{t)  =  ex.p{eV't)S^ .  (26) 

Finally,  we  obtain  the  approximate  solution  as 

p(i)  =  v(l)p(0).  (27) 

As  an  illustration,  in  Figure  1  we  show  components  pi{t) 
and  P2{t)  of  the  solution  above  for  the  initial  condition 
Pi{0)  =  S2i,  and  e  =  0.01.  We  can  see  that  after  the  tran¬ 
sient  time  (17)  has  elapsed  we  obtain  a  good  agreement 
between  the  exact  and  the  approximate  solution  to  this 
example  problem. 

To  further  support  our  results,  we  randomly  generate  a 
large  number  of  master  equations  with  similar  near  block- 
diagonal  structure  and  compare  their  exact  solutions  to 
the  approximate  solutions  obtained  using  our  approach. 
The  numerical  results  presented  in  Figure  2  show  that  the 
approximation  error  is  controlled  by  the  small  parameter 
e. 


IV.  PAP  SWITCH 

Pili  are  small  hair-like  structures  that  enable  bacte¬ 
ria  to  bind  to  epithelial  cells  and  thereby  significantly 


The  Pap  Operon 


Dam  methylation  (GATC)  sites 

n  ryn  n 

U  U  U  U  LI  U 

4  5  6  1  2  3 

Lrp  Binding  Sites 


The  Chemical  Regulators 


Lrp-Leucine  responsive 
regulatory  protein 
Dam-DNA  adenine 
methylase 


Diagram  of  an  ON  State 

•  Lrp  is  bound  at  sites  4-5-6 

•  Dam  has  methylated  site  2 


FIG.  3:  Schematic  of  the  pap  operon  (top),  key  regulatory 
components  of  the  Pap  switch  (middle)  and  diagram  of  the 
operon  in  its  on  state  (bottom). 


increase  the  bacteria’s  ability  to  infect  host  organisms. 
However,  pili  expression  comes  at  cost  to  the  bacteria,  as 
the  production  of  pili  requires  a  large  portion  of  the  cel¬ 
lular  energy.  Whether  or  not  E.  coli  are  piliated  depends 
upon  the  regulation  of  genes  such  as  the  pyelonephritis- 
associated  pili  (pap)  genes.  Here  we  study  a  simplified 
version  of  the  PAP  switch  modeF^,  which  analyzes  the 
regulatory  network  responsible  for  controlling  one  type 
of  pili. 

Recent  experiments^^’^^  have  identified  two  transcrip¬ 
tion  factors  that  affect  expression  of  the  pap  gene,  and  six 
binding  sites  for  the  two.  The  transcription  factors  are 
DNA  adenine  methylase  (Dam)  and  leucine  responsive 
regulatory  protein  (Lrp).  Dam  binds  and  applies  methyl 
groups  to  GATC  sites  at  2  and  5,  as  shown  in  Figure 
3.  This  Dam  methylation  is  an  irreversible  and  relatively 
slow  process.  On  the  other  hand,  Lrp  binds  cooperatively 
to  three  adjacent  sites  at  a  time,  either  1-2-3  or  4-5-6 
(Figure  3).  These  reactions  are  fast  and  reversible.  Lrp 
binding  also  inhibits  Dam  methylation.  Altogether  this 
makes  16  possible  states  in  which  the  PAP  switch  can  be 
found.  These  are  schematically  described  by  the  network 
model  shown  in  Figure  4.  In  this  model  pap  transcription 
occurs  only  in  the  state  11  when  Lrp  is  bound  to  sites 
1-2-3  and  site  5  is  methylated.  We  assume  that  cell  repli¬ 
cation  always  resets  the  system  to  state  1.  A  solution  to 
the  chemical  master  equation  for  this  problem  gives  the 
time  evolution  of  the  probability  of  finding  the  system 
in  each  state  including  state  11,  which  is  proportional  to 
the  probability  of  pap  gene  expression. 

Since  the  two  transcription  factors  bind  at  significantly 
different  rates,  a  smart  bookkeeping  practice  would  be 
to  record  Lrp  binding  propensities  in  the  matrix  H  and 
methylation  propensities  in  V  as  defined  in  (8).  With  a 
convenient  labeling  scheme,  as  shown  in  Figure  4,  we  can 
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FIG.  4:  PAP  switch  schematic  diagram. 


express  in  a  simple  block  diagonal  form, 

/JJi  0  0  0  \ 

0  172  0  0 

0  0  iJa  0  ■ 

\  0  0  0  Hi  / 


(28) 


Recent  experimental  data^^  reveals  that  the  propensi- 
ties  of  Lrp  binding  at  sites  4-5-6  depend  strongly  on  the 
methylation  pattern  of  site  5,  while  propensities  of  Lrp  at 
sites  1-2-3  do  not  significantly  depend  upon  the  methy- 
lation  pattern  of  site  2.  Thus  we  find  that  there  are  only 
two  distinct  blocks  as 


Hi  =Hs  = 


and 


Ho  =  Ha  = 


-9500 

6.8 

0.09 

0 

9270 

-18.4 

0 

0.09 

230 

0 

-463.29 

6.79 

0 

11.6 

463.2 

-6.88 

-9500 

62 

0.09 

0  ' 

9270 

-73 

0 

0.09 

230 

0  - 

-463.29 

61.76 

0 

11 

463.2  - 

-61.85  . 

(29) 


(30) 


Leading  eigenvalues  for  both,  Hi  and  H2,  are  zero,  while 
the  next  largest  eigenvalue  is  of  order  A5  ^  —10.  On  the 
other  hand  we  estimate  that  all  methylation  propensities 
have  the  same  value  e  =  0.17.  Following  our  labeling 
scheme  (Figure  4)  the  nonzero  entries  of  matrix  V  are 
then:  Vi^i  =  —2,  ¥2,2  =  V3,3  =  ^5,5  =  Vjj  =  Vg^g  = 
^10,10  =  — and  Vsp  =  Vgp  =  V6,2  =  = 

hi3,g  =  Vi4po  =  Vi5j  =  1-  Therefore,  all  we  need  to  con¬ 
struct  the  matrix  are  the  right  eigenvectors  vi  and 
V2  that  correspond  to  the  zero  eigenvalues  of  Hi  and  iL2, 
respectively.  Following  the  footsteps  outlined  in  Section 
III  we  reduce  the  PAP  switch  model  to  a  4-dimensional 
system  and  carry  out  calculation  for  the  probability  pn, 
which  is  proportional  to  the  PAP  transcription  probabil¬ 
ity. 

The  PAP  switch  model  we  presented  here  is  simple 
enough  to  be  integrated  directly  so  we  can  compare  re¬ 
sults  for  the  full  system  and  the  reduced  system.  As  we 


FIG.  5:  Time  evolution  of  pap  gene  expression  probability. 
Initially  no  transcription  factors  are  bound  the  pap  operon,  so 
the  initial  condition  is  pi(0)  =  1  and  Pi^i  =  0.  The  transient 
time  17  is  less  than  1  in  our  time  units. 


show  in  the  Figure  5,  all  the  important  information  about 
the  system’s  behavior  is  preserved  in  the  reduced  model. 

This  model  predicts  a  short  time  lag  between  repli¬ 
cation  and  Pap  production,  since  methylation  of  site  2 
must  occur  before  pap  expression.  Further,  since  Dam 
methylation  at  5  prohibits  expression,  if  the  cell  waits 
too  long  to  decide  to  switch  “on” ,  it  will  most  probably 
miss  its  chance  and  remain  “off” .  Thus,  a  newly  created 
E.  coli  cell  will  most  likely  express  the  pap  gene  at  some 
point  shortly  after  replication.  Probability  of  expressing 
pili  drops  significantly  at  later  times  and  cell  resources 
are  used  for  other  functions,  such  as  initiating  the  next 
replication  cycle. 


V.  HEAT  SHOCK  SYSTEM 

Through  evolution  all  living  organisms  have  developed 
mechanisms  for  dealing  with  environmental  stress.  One 
such  mechanism  is  cellular  heat  shock  response.  In¬ 
creased  temperature  causes  proteins  inside  the  cell  to 
misfold  and  thereby  loose  their  functionality.  In  re¬ 
sponse,  the  cell  produces  heat  shock  proteins,  most  no- 
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FIG.  6:  (a)  Two  dimensional  lattice  representing  possible 

states  and  transitions  in  the  heat  shock  model.  Here  si  and 
S3  are  populations  of  cr32-DnaK  and  cr32-RNAP  compounds, 
respectively,  while  S2  is  the  population  of  free  (732  molecules. 
Reactions  si  ^  S2,  are  represented  by  bidirectional  horizontal 
arrows  and  reactions  S2  ^  S3  is  represented  with  vertical 
arrows.  The  total  number  of  732  is  constant  (in  this  example 
Si  +  S2  +  S3  =  5),  so  the  chemical  state  of  the  system  is 
uniquely  dehned  by  si  and  S3  alone,  (b)  The  same  lattice 
after  applying  the  Finite  State  Projection.  Unlikely  states 
have  been  aggregated  into  a  single  sink  state. 


tably  molecular  chaperones  and  proteases.  Molecular 
chaperones  help  refold  deformed  proteins  and  restore 
their  original  function.  Proteases,  on  the  other  hand, 
help  degrade  damaged  proteins,  hopefully  before  they 
trigger  unwanted  chemical  reactions  in  the  cell.  At  the 
crux  of  the  heat  shock  response  mechanism  in  E.  coli  is 
the  synthesis  of  cr32-RNAP  complex^^.  Here  we  study  a 
simplified  model  for  cr32-RNAP  synthesis  using  our  time 
scale  separation  method  in  conjunction  with  a  Finite 
State  Projection.  The  analysis  of  the  entire  heat  shock 
response  mechanism  in  E.  coli  is  rather  complicated,  and 
is  beyond  the  scope  of  this  paper. 

At  normal  physiological  temperatures  (732  protein  is 
found  almost  exclusively  in  a  complex  cr32-DnaK.  As  the 
temperature  increases  this  complex  becomes  less  stable 


FIG.  7:  Probability  that  no  cr32-RNAP  molecule  has  been 
sythesized  in  the  heat  shock  toy  model. 


and  there  is  a  non-negligible  probability  of  finding  free 
<732  inside  the  cell.  The  free  <732  then  can  combine  with 
RNA  polymerase  through  what  can  be  considered  an  irre¬ 
versible  reaction  to  form  a  (732-RNAP  complex.  In  turn, 
cr32-RNAP  initiates  transcription  of  genes  encoding  heat 
shock  proteins.  This  regulatory  mechanism  can  be  sum¬ 
marized  in  a  simple  set  of  reactions, 


Si  ^  S2  ^  S3,  (31) 

where  5i,  S2  and  S3  correspond  to  the  cr32-DnaK  complex, 
the  <732  heat  shock  regulator  and  the  <732-RNAP  complex, 
respectively.  This  model  of  the  heat  shock  subsystem 
has  been  analyzed  before  using  various  computational 
methods  including  Monte  Carlo  implementations'’^^. 

In  the  biological  system,  the  relative  rates  of  the  reac¬ 
tions  are  such  that  the  reaction  from  S2  to  Si  is  by  far 
the  fastest,  and  (732  molecules  infrequently  escape  from 
DnaK  long  enough  to  form  the  cr32-RNAP  complex.  The 
purpose  of  this  mechanism  is  to  strike  a  balance  between 
fixing  the  damage  produced  by  heat  and  saving  the  cell’s 
resources,  as  a  significant  portion  of  cell  energy  is  con¬ 
sumed  when  producing  heat  shock  proteins.  The  optimal 
response  to  the  heat  shock  is  not  massive,  but  measured 
production  of  heat  shock  proteins,  which  leaves  sufficient 
resources  for  other  cellular  functions.  We  use  the  follow¬ 
ing  set  of  parameters  values  for  the  reaction  rates^’^^: 

Cl  =  10,  C2  =  4  X  10"^,  C3  =  2.  (32) 

For  simplicity,  in  our  model  we  assume  that  the  total 
number  of  <732  -  free  or  in  compounds  -  is  constant,  so 
that  Si  +  52  +  53  =  const.  With  this  constraint  the  reach¬ 
able  states  of  this  three  species  problem  can  be  repre¬ 
sented  on  a  two  dimensional  lattice. 

For  illustrative  purposes.  Figure  6a,  shows  one  such 
lattice  for  an  initial  condition  of  si  =  5  and  S2  =  S3  =  0. 
Here,  the  total  population  is  fixed  at  five,  and  there  is  a 
total  of  21  reachable  states. 

We  first  apply  the  Finite  State  Projection.  We  esti¬ 
mate  that  all  states  where  S2  >  2  or  S3  >  2  are  unlikely 
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to  be  reached  in  a  short  time,  so  we  aggregate  them  into 
a  sink  node  as  shown  in  Figure  6b,  thereby  reducing  this 
to  a  ten  state  problem.  From  the  transitions  to  the  ag¬ 
gregated  state,  we  find  a  strict  upper  bound  on  the  error 
introduced  by  such  an  approximation.  For  our  set  of  pa¬ 
rameters  the  approximation  the  error  is  guaranteed  to  be 
below  8%  for  any  time  t  ^  500. 

Next,  we  further  reduce  this  system  by  applying  time 
scale  separation.  Elements  of  the  matrix  ApsPi  which 
defines  the  master  equation  for  the  system  obtained  after 
the  Finite  State  Projection,  can  be  read  off  of  Figure  6b. 
In  accordance  with  our  bookkeeping  practice  we  write 
Apsp  =  H  and  record  all  reversible  reactions  si  ^ 
52  in  the  matrix  and  all  other  reactions,  including 
S2  ss  and  transitions  to  the  aggregated  state,  in  the 
matrix  V.  By  doing  so  we  ensure  that  all  fast  reactions 
are  contained  in  H.  Note  that  there  is  no  unique  way  to 
separate  fast  and  slow  reactions  and  we  chose  this  one  for 


its  simplicity.  Matrix  H  has  a  block  diagonal  structure 
/Hs 


H  = 


Ho 


Hi 


(33) 


where  each  block 


-{k  +  2)ci  C2  0 

Hk  =  (  (fc  +  2)ci  ~{k  +  l)ci  —  C2  2c2 
0  {k  +  l)ci  -2c2 


(34) 


corresponds  to  a  row  of  states  in  Figure  6b.  The  zero  in 
the  last  row  is  just  a  scalar,  and  it  corresponds  to  the  ag¬ 
gregated  state.  The  matrix  eV  is  made  up  of  irreversible 
reactions  (vertical  transitions  in  Figure  6b)  and  therefore 
has  a  lower  triangular  form: 
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For  the  reaction  rates  above,  the  first  four  eigenvalues 
of  H  are  zero,  and  the  rest  have  negative  real  parts  each 
with  magnitude  of  order  10^  or  larger,  suggesting  that  the 
truncation  in  (16)  is  indeed  valid  for  this  problem.  There¬ 
fore  the  dynamics  of  9-dimensional  system  obtained  by 
the  Finite  State  Projection  can  be  well  approximated  by 
a  system  of  only  four  linear  ordinary  differential  equa¬ 
tions.  By  applying  algorithm  from  Section  III  A  we  find 
that  the  time  scale  separation  introduces  error  of  order 
10“^  with  respect  to  the  solution  obtained  by  Finite  State 
Projection  alone.  The  transient  time  (17)  is  estimated  to 
be  2  X  10“"^,  and  is  negligible  considering  the  time  interval 
of  interest. 

In  Figure  7  we  compare  results  obtained  by  solving  the 
full  system  directly,  using  Finite  State  Projection  alone 
and  using  Finite  State  Projection  and  time  scale  sepa¬ 
ration  combined.  The  figure  shows  how  probability  of 
having  no  cr32-RNAP  complex  in  the  cell  decreases  with 
time.  All  three  results  are  in  a  good  agreement  as  our 
calculations  predicted.  Note  that  the  Finite  State  Pro¬ 
jection  error  in  this  example  is  much  smaller  than  the 
estimated  upper  bound. 

The  advantage  of  combining  the  Finite  State  Projec¬ 


tion  and  time  scale  separation  becomes  obvious  if  we  con¬ 
sider  a  more  realistic  and  much  larger  problem  with  the 
same  reactions  but  with  initial  conditions  si  =  2000  and 
<^2  =  <^3  =  0.  In  this  case  there  are  2,001,000  reach¬ 
able  states,  and  the  full  chemical  master  equation  is  too 
large  to  be  tackled  directly.  However,  by  applying  the  Fi¬ 
nite  State  Projection  we  find  that  truncating  every  state 
where  53  ^  340  and  52  ^  15  introduces  1-norm  error  that 
is  less  than  10“^  for  times  t  ^  3005.  The  resulting  matrix 
A  is  of  size  5457  x  5457,  and  has  near  block  diagonal  form 
(8)  similar  to  the  example  in  Figure  6.  Its  block  diago¬ 
nal  part  H  contains  341  irreducible  blocks  each  with  16 
rows  and  columns.  Same  as  in  the  previous  example  the 
leading  nonzero  eigenvalue  of  H  has  a  negative  real  part 
of  magnitude  10^,  so  the  system  can  be  reduced  to  a  342 
state  model  using  the  time  scale  separation  algorithm. 
Should  we  apply  time  scale  separation  directly  to  the  full 
system,  not  only  the  amount  of  the  computation  would 
be  significantly  larger,  but  there  would  be  no  simple  way 
to  obtain  the  near  block  diagonal  structure  as  it  was  the 
case  in  the  previous  example. 

The  solution  to  this  problem  shows  how  the  number 
of  compounds  cr32-RNAP  grows  in  time  if  the  temper- 
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FIG.  8:  Probability  distribution  for  S3  calculated  at  three 
different  times.  The  truncated  solution  (dots)  approximates 
well  the  solution  to  the  full  system  (solid  lines). 

ature  is  constant  and  slightly  above  normal  physiologi¬ 
cal  level.  This  number  is  proportional  to  the  number  of 
heat  shock  proteins  produced  in  the  cell.  With  the  Finite 
State  Projection  solution,  we  have  computed  the  proba¬ 
bility  distribution  for  S3  at  three  times  t  =  100,  200,  and 
300.  Figure  8  illustrates  the  probability  distribution  of 
species  cr32-RNAP  (53)  at  these  instances  in  time.  The 
Finite  State  Projection  solution,  shown  with  solid  lines, 
takes  approximately  220  seconds  to  compute.  The  dots 
represent  the  probability  distribution  of  S3  as  computed 
using  our  time  scale  separation  method  applied  atop  of 
the  Finite  State  Projection.  The  difference  between  the 
two  results  is  indistinguishable,  but  the  total  computa¬ 
tional  effort  for  the  reduced,  342  state  model  is  200-fold 
less;  the  truncated  model  takes  only  about  a  second  to 
solve  the  problem. 

VI.  CONCLUSION 

Until  recently,  it  was  thought  that  the  chemical  master 
equation  could  not  be  analytically  solved  except  for  the 
most  trivial  systems.  Previous  work  on  the  Finite  State 
Projection  demonstrated  that  for  many  biological  sys¬ 
tems,  bulk  system  reductions  could  bring  models  closer 
into  the  fold  of  solvable  problems.  Here  we  have  shown 
that  the  Finite  State  Projection  method  can  be  further 
enhanced  when  solving  the  chemical  master  equation  for 
systems  involving  multiple  time  scales.  In  combination 
with  Finite  State  Projection  method,  we  have  shown  that 
our  algorithm,  based  upon  singular  perturbation  theory, 
provides  a  powerful  computational  tool  for  studying  in- 
tracelular  processes  and  gene  regulatory  networks. 

Similar  problems  were  studied  earlier  with  specially 
designed  Monte  Carlo  implementations^’^  or  hybrid 
methods^.  In  contrast  to  these,  our  method  does  not 
require  random  number  generation,  and  its  accuracy  is 
given  a  priori.  A  further  advantage  of  our  method  is  its 
ease  of  implementation  and  the  speed  of  computations. 


The  proposed  algorithm  is  particularly  fast  when  imple¬ 
mented  on  systems  for  which  there  are  strict  means  of 
separating  slow  and  fast  reactions. 

The  Finite  State  Projection  and  our  time  scale  sepa¬ 
ration  approach  also  provide  valuable  insight  as  to  how 
one  may  further  deal  with  the  bewildering  complexity 
that  intracellular  processes  exhibit.  First,  cellular  pro¬ 
cesses  are  limited  by  cell  size  and  available  energy.  It 
is  then  plausible  that  the  main  features  of  intracellular 
dynamics  can  be  captured  in  a  relatively  small  subset  of 
the  state  space,  as  the  results  obtained  by  Finite  State 
Projection  suggest.  Another  typical  feature  of  intracellu¬ 
lar  processes  is  that  they  are  composed  of  reactions  that 
take  place  on  different  timescales.  Depending  on  the  ob¬ 
servation  time  of  interest,  some  of  these  reactions  can  be 
neglected,  while  some  will  contribute  only  through  their 
averages.  Preliminary  success  with  our  approach  gives 
us  a  hope  that  relatively  simple  models  for  intracellular 
processes  can  be  tailored  when  a  region  in  the  state  space 
and  observation  time  of  interest  are  known. 

Of  course,  one  can  easily  envision  that  additional 
model  reductions  may  be  possible  to  even  further  en¬ 
hance  the  power  of  both  the  Finite  State  Projection  and 
the  time  scale  separation  approach.  Indeed  some  reduc¬ 
tions  based  upon  control  theory^^  are  already  becoming 
apparent.  Also,  in  our  computations  we  have  used  off  the 
shelf  numerical  routines  for  eigensystem  calculations  and 
matrix  exponentiation.  Further  improvements  in  compu¬ 
tational  speed  can  be  achieved  if  these  routines  are  op¬ 
timized  for  matrices  which  define  master  equations  and 
their  special  properties.  We  intend  to  investigate  these 
possibilities  in  the  future. 
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APPENDIX 

Singular  perturbation  theory  has  been  extensively 
studied  in  various  literature.  However,  most  of  the  litera¬ 
ture  in  this  area  is  of  wide  scope  and  often  very  technical. 
In  order  to  spare  the  reader  some  time,  here  we  present  a 
heuristic  argument,  which  provides  a  mathematical  jus¬ 
tification  for  our  method,  while  keeping  technicalities  at 
minimum.  For  rigorous  proofs  interested  reader  may 
want  to  consult  for  example^"^’^^. 

Consider  a  weakly  perturbed  linear  A-dimensional  sys- 
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tern  described  by  (11) 


N 

Xi  =  X^Xi  -|-  6  ^  ^  ^ij^j  •) 

i=i 


interval  of  interest.  These  equations  are  linear  and  hence 
can  be  solved  analytically,  but  let  us  take  an  extra  step 
here  and  expand  the  solution  to  (A.l)  in  powers  of  e 

Xi{t)  =  xf\t)  +  +  e^xf\t)  +  . . .  (A. 6) 


where  =  0  for  i  ^  m,  and  has  negative  real  part  for 
i  >  m.  We  want  to  find  a  near  identity  coordinate  trans¬ 
formation  (12)  that  would  remove  as  many  0(e)  terms 
as  possible  from  (A.l)  and  “push”  them  to  higher  orders 
in  e.  After  substituting  x  =  (/  +  eG)y  in  (A.l)  we  get 

N  N 

iji  —  XiHi  +  e  ^  ^  ^ijVj  ^  ^  ^ 
i=i  i=i 

N 

-^^Xi  GijUj  +  O(e^) 
i=i 


By  substituting  this  expression  into  (A.l)  and  grouping 
same  orders  in  e  we  get  series  of  equations 

e^  :  xf\t)  =  XiX^P  {t) 

N 

=  Xix\^\t)  +  Vijxf\t) 

N 

:  xf\t)  =  \ixf\t)  +  y]  Vijxf\t) 
i=i 

which  we  can  solve  in  a  straightforward  way  to  obtain 


By  equating  all  0(e)  terms  to  zero  we  find 

N 

(Tdi  -  GijXj  +  XiGij)  yj  =  0, 

i=i 


N  .t 

V  Vijxf\Q)  / 

i=i 


and  by  solving  for  we  obtain 


Gi,= 


Vi, 


\j  Aj 


(A.4) 


Therefore,  we  can  always  find  Gij  except  when  Xi  =  Xj. 
In  other  words,  all  nonresonant  terms  can  be  removed 
through  0(e)  from  (A.l)  by  a  near  identity  transforma¬ 
tion  (12).  In  our  method  we  are  interested  in  separating 
slow  and  fast  processes  in  the  system,  so  we  shall  define 
matrix  G  in  (12)  as  follows: 


Gij  - 


y^J 

Xj  A^ 
0 


i  ^  m  <  j 
otherwise 


(A.5) 


By  substituting  this  expression  for  G  in  (A. 2)  we  find 
that 


m 

y*  =  e  y]  Vijyj  +  O(e^)  i  ^  m 

i=i 

N 

m  =  XiVi  +  e  y]  VijVj  +  O(e^)  i  >  m 
i=i 

We  observe  that  first  m  equations  decouple  from  the  rest 
of  the  system  through  0(e),  and  can  be  solved  indepen¬ 
dently  after  truncating  higher  order  terms.  Furthermore, 
the  near  identity  transformation  (12)  does  not  introduce 
any  new  0{e)  terms  to  the  first  m  equations,  so  it  is  es¬ 
sentially  just  a  truncation  of  all  Xi^rn  terms  from  (A.l). 
We  do  not  need  to  calculate  G  and  perform  transforma¬ 
tion  (12)  as  such  transformation  is  guaranteed  to  exist. 

It  remains  to  show  that  the  solution  to  truncated  sys¬ 
tem  (13)  will  be  0(e)  close  to  solution  to  (A.l)  on  a  time 


Let  us  first  consider  equations  i  ^  m.  The  solution  to 
(A.l)  through  0(e)  then  can  be  written  as 

m 

Xi{t)  =  Xi{0)  +  ey^VijXj{0)t 

i=i 

pAjt  _  1 

j=m-\-l  ^ 

where  we  substituted  ^^(O)  =  +  ex^^^^O)  +  O(e^) 

and  A^  =  0.  Since  the  system  has  one  stable  steady  state 
solution  the  series  above  must  converge  for  all  times.  The 
first  two  terms  in  the  expansion  above  are  equal  to  the 
first  two  terms  in  the  expansion  of  the  solution  (14)  for 
the  truncated  system.  Therefore,  for  yi{0)  =  ^^(0)  it  is 


\xi{t)  -yi{t)\ 


N 

E 

j=m-\-l 


^Xjt 


-V-  -x^^^ 

VljXj 


(0) 


+  0(6^). 


In  the  expression  above  all  Xj  <  0,  therefore  \xi{t)  — 
yi{t)\  =  0(e)  holds  for  all  t  >  0.  Since  the  expression  for 
Xi{t)  is  convergent  series  and  ^^(0)  are  linearly  indepen¬ 
dent,  we  conclude  that  yi{t)  must  also  have  a  fixed  point 
solution,  which  is  0(e)  close  to  the  solution  of  the  full 
system. 

Next  we  consider  equations  in  (A.l)  where  i  >  m.  The 
solution  to  these  can  be  expanded  in  terms  of  e  as 

^Xit  ^  ^ 

Xi{t)  =  e^‘*a;i(0) +  e^— - 

*  i=i 

N 

^  ^  VijXjifi)  -T  O(e^) 
i=m+i 
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Our  truncation  algorithm  (Sec.  Ill  A)  sets  all  yi{t)  = 
0,  so  initially  the  difference  between  full  and  truncated 
solution  is  whatever  the  initial  condition  x^(0)  is,  and  it 
can  be  larger  than  0(e).  However,  in  the  limit  case 


lim  \xi{t)\ 


^  m 


i=i 


0(e2)  (A.7) 


That  means  the  truncation  introduces  0(e)  error  to  the 
asymptotic  solution.  Larger  errors  may  occur  only  during 
the  finite  transient  time  0  <  t  <  T(e),  where  T(e)  is  given 
in  (17).  One  can  verify  this  by  substituting  right  hand 
side  of  (17)  for  time  in  the  solution  x^>^(t)  above. 
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We  shall  assume  through  the  rest  of  the  paper  that  N 
is  finite.  Eor  systems  in  which  N  is  infinite,  the  Einite 
State  Projection  could  be  applied  to  find  a  finite  system, 
which  approximates  the  original  to  an  arbitrary  degree  of 
accuracy. 


